Load libraries need to process the data

# load library
library(raster)
## Loading required package: sp
library(rhdf5)
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3
# check working directory
getwd()
## [1] "/Users/cassondrawalker/Documents/data/NEONDI_2016/monday_pm"

Open data file

# Open file object by defining file path 
# "tab" key allows you to select an element in a folder
# ../up one folder

f <- "../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

# View h5 structure

h5ls(f)
##    group                                 name       otype  dclass
## 0      /                     ATCOR_Input_File H5I_DATASET  STRING
## 1      /                 ATCOR_Processing_Log H5I_DATASET  STRING
## 2      /                Aerosol Optical Depth H5I_DATASET INTEGER
## 3      /                               Aspect H5I_DATASET   FLOAT
## 4      /                          Cast Shadow H5I_DATASET INTEGER
## 5      / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6      /                 Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7      /                  Illumination Factor H5I_DATASET INTEGER
## 8      /                          Path Length H5I_DATASET   FLOAT
## 9      /                   Processing Version H5I_DATASET  STRING
## 10     /                          Reflectance H5I_DATASET INTEGER
## 11     /                Shadow_Processing_Log H5I_DATASET  STRING
## 12     /                      Sky View Factor H5I_DATASET INTEGER
## 13     /               Skyview_Processing_Log H5I_DATASET  STRING
## 14     /                                Slope H5I_DATASET   FLOAT
## 15     /          Slope_Aspect_Processing_Log H5I_DATASET  STRING
## 16     /                  Solar Azimuth Angle H5I_DATASET   FLOAT
## 17     /                   Solar Zenith Angle H5I_DATASET   FLOAT
## 18     /                    Surface Elevation H5I_DATASET   FLOAT
## 19     /                 Visibility Index Map H5I_DATASET INTEGER
## 20     /                   Water Vapor Column H5I_DATASET INTEGER
## 21     /             coordinate system string H5I_DATASET  STRING
## 22     /                       flightAltitude H5I_DATASET   FLOAT
## 23     /                        flightHeading H5I_DATASET   FLOAT
## 24     /                           flightTime H5I_DATASET   FLOAT
## 25     /                                 fwhm H5I_DATASET   FLOAT
## 26     /                             map info H5I_DATASET  STRING
## 27     /              to-sensor azimuth angle H5I_DATASET   FLOAT
## 28     /               to-sensor zenith angle H5I_DATASET   FLOAT
## 29     /                           wavelength H5I_DATASET   FLOAT
##                dim
## 0                1
## 1                1
## 2    544 x 578 x 1
## 3    544 x 578 x 1
## 4    544 x 578 x 1
## 5    544 x 578 x 1
## 6    544 x 578 x 1
## 7    544 x 578 x 1
## 8    544 x 578 x 1
## 9                1
## 10 544 x 578 x 426
## 11               1
## 12   544 x 578 x 1
## 13               1
## 14   544 x 578 x 1
## 15               1
## 16           1 x 1
## 17           1 x 1
## 18   544 x 578 x 1
## 19   544 x 578 x 1
## 20   544 x 578 x 1
## 21               1
## 22         5732053
## 23         5732053
## 24         5732053
## 25         426 x 1
## 26               1
## 27   544 x 578 x 1
## 28   544 x 578 x 1
## 29         426 x 1

Import Spatial Information

# Import spatial file

mapInfo <- h5read(file = f,
                name = "map info",
                read.attributes = TRUE)
mapInfo
## [1] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## attr(,"Description")
## [1] "Basic Map information for envi style programs"
# read attributes brings in the metadata as well so you can access it

Import Special Attribute Values - Reflectance Metadata

# read reflectance data attributes

reflInfo <- h5readAttributes(file = f,
                           name = "Reflectance")
reflInfo
## $DIMENSION_LABELS
## [1] "Wavelength" "Line"       "Sample"    
## 
## $Description
## [1] "Atmospherically corrected reflectance."
## 
## $`Scale Factor`
## [1] 10000
## 
## $Unit
## [1] "unitless. Valid range 0-1."
## 
## $`data ignore value`
## [1] "15000.0"
# Define Scale Factor
scaleFactor <- reflInfo$`Scale Factor`

# Define No Data Value

noDataValue <- reflInfo$`data ignore value`

# Look at strucutre of data
str(scaleFactor)
##  num 10000
str(noDataValue)
##  chr "15000.0"
# Redefine no data value as numeric, not character
noDataValue <- as.numeric(reflInfo$`data ignore value`)

Import Data Dimensions

# open file for viewing-connection to file made (f for file)
fid <- H5Fopen(f)
fid
## HDF5 FILE
##         name /
##     filename 
## 
##                                    name       otype  dclass
## 0  ATCOR_Input_File                     H5I_DATASET STRING 
## 1  ATCOR_Processing_Log                 H5I_DATASET STRING 
## 2  Aerosol Optical Depth                H5I_DATASET INTEGER
## 3  Aspect                               H5I_DATASET FLOAT  
## 4  Cast Shadow                          H5I_DATASET INTEGER
## 5  Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6  Haze-Cloud-Water Map                 H5I_DATASET INTEGER
## 7  Illumination Factor                  H5I_DATASET INTEGER
## 8  Path Length                          H5I_DATASET FLOAT  
## 9  Processing Version                   H5I_DATASET STRING 
## 10 Reflectance                          H5I_DATASET INTEGER
## 11 Shadow_Processing_Log                H5I_DATASET STRING 
## 12 Sky View Factor                      H5I_DATASET INTEGER
## 13 Skyview_Processing_Log               H5I_DATASET STRING 
## 14 Slope                                H5I_DATASET FLOAT  
## 15 Slope_Aspect_Processing_Log          H5I_DATASET STRING 
## 16 Solar Azimuth Angle                  H5I_DATASET FLOAT  
## 17 Solar Zenith Angle                   H5I_DATASET FLOAT  
## 18 Surface Elevation                    H5I_DATASET FLOAT  
## 19 Visibility Index Map                 H5I_DATASET INTEGER
## 20 Water Vapor Column                   H5I_DATASET INTEGER
## 21 coordinate system string             H5I_DATASET STRING 
## 22 flightAltitude                       H5I_DATASET FLOAT  
## 23 flightHeading                        H5I_DATASET FLOAT  
## 24 flightTime                           H5I_DATASET FLOAT  
## 25 fwhm                                 H5I_DATASET FLOAT  
## 26 map info                             H5I_DATASET STRING 
## 27 to-sensor azimuth angle              H5I_DATASET FLOAT  
## 28 to-sensor zenith angle               H5I_DATASET FLOAT  
## 29 wavelength                           H5I_DATASET FLOAT  
##                dim
## 0  1              
## 1  1              
## 2  544 x 578 x 1  
## 3  544 x 578 x 1  
## 4  544 x 578 x 1  
## 5  544 x 578 x 1  
## 6  544 x 578 x 1  
## 7  544 x 578 x 1  
## 8  544 x 578 x 1  
## 9  1              
## 10 544 x 578 x 426
## 11 1              
## 12 544 x 578 x 1  
## 13 1              
## 14 544 x 578 x 1  
## 15 1              
## 16 1 x 1          
## 17 1 x 1          
## 18 544 x 578 x 1  
## 19 544 x 578 x 1  
## 20 544 x 578 x 1  
## 21 1              
## 22 5732053        
## 23 5732053        
## 24 5732053        
## 25 426 x 1        
## 26 1              
## 27 544 x 578 x 1  
## 28 544 x 578 x 1  
## 29 426 x 1
# open the reflectance dataset-make connection to data (d for data)
did <- H5Dopen(fid,"Reflectance")
did
## HDF5 DATASET
##         name /Reflectance
##     filename 
##         type H5T_STD_I16LE
##         rank 3
##         size 544 x 578 x 426
##      maxsize 544 x 578 x 426
# data is read as column,row,bands which is different from HDF5 Viewer

# grab dataset dimensions (x,y,z data) s for structure
sid <- H5Dget_space(did)
sid
## HDF5 DATASPACE
##         rank 3
##         size 544 x 578 x 426
##      maxsize 544 x 578 x 426
# Import the dimesions of the data as column, row, band
dims <- H5Sget_simple_extent_dims(sid)$size
dims
## [1] 544 578 426
# close all open connections-otherwise could overwrite data
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)

Read in Reflectance Data

# Extract slice of H5 file
# Import BAND 56 index is list of things to import (column,rows,bands)
 b56 <- h5read(file = f,
               name = "Reflectance",
               index = list(1:dims[1],1:dims[2],56))

Convert Data to Matrix

# convert to matrix
b56 <- b56[,,1]

# plot data
image(b56)

image(log(b56),
      main = "log transformed data")

# Look at distribution of data
hist(b56)

Clean the data

# assign no data value to object-will not be used in computations
b56[b56 == noDataValue] <- NA

# apply scale factor
b56 <- b56/scaleFactor

hist(b56)

Transpose Data (Rows and Columns)

# transpose the row and colums
b56 <- t(b56)

image(log(b56))

Create Spatial Extent

# Find data structure
str(mapInfo)
##  chr [1(1d)] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
##  - attr(*, "Description")= chr "Basic Map information for envi style programs"
class(mapInfo)
## [1] "array"
# split out Map Info object
mapInfo <- strsplit(mapInfo, ",")
# Take out heirarchy strucutre in the data
mapInfo <- unlist(mapInfo)

# look at data structure
str(mapInfo)
##  chr [1:11] "UTM" "1" "1" "325963.0" "4103482.0" ...
xMin <- as.numeric(mapInfo[4])
yMin <- as.numeric(mapInfo[5])

# to find the r coordinate for  x and y max  take min+(dimensions*resolution)
# to get spatial resolution
xres <- as.numeric(mapInfo[6])
yres <- as.numeric(mapInfo[7])

# bring in the coordinates to convert r coordinates to spatial coordinates
xMax <- xMin + (dims[1] * xres)
yMax <- yMin + (dims[2] * yres)

Create Extent Object

# Create extent (order xmin, xmax, ymin, ymax) 
rasExt <- extent(xMin, xMax, yMin, yMax)
rasExt
## class       : Extent 
## xmin        : 325963 
## xmax        : 326507 
## ymin        : 4103482 
## ymax        : 4104060

Create actual raster object

# raster with coordinate system extent
b56r <- raster(b56,
               crs=CRS("+init=epsg:32611"))

# assign extent to raster
extent(b56r) <- rasExt
b56r
## class       : RasterLayer 
## dimensions  : 578, 544, 314432  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326507, 4103482, 4104060  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0033, 0.5678  (min, max)
plot(b56r)

Import NEON functions

# install devtools
# install.packages("devtools")
library(devtools)

install_github("lwasser/neon-aop-package/neonAOP")
## Skipping install for github remote, the SHA1 (5b9a37fe) has not changed since last install.
##   Use `force = TRUE` to force installation
library(neonAOP)
# use open band function to view data
b55 <- open_band(f,
                 bandNum = 55,
                 epsg = 32611)
b55 
## class       : RasterLayer 
## dimensions  : 578, 544, 314432  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326507, 4102904, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0029, 0.5622  (min, max)
# plot data
plot(b55)

# import several bands (R band,G band,B band)
bands <-c(58, 34, 19)
epsg <- 32611

# create raster stack
RGBstack <- create_stack(f,
                         bands = bands,
                         epsg = epsg)

RGBstack
## class       : RasterStack 
## dimensions  : 578, 544, 314432, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326507, 4102904, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :    X58,    X34,    X19 
## min values  : 0.0030, 0.0027, 0.0009 
## max values  : 0.5680, 0.4761, 0.3808
plot(RGBstack)

plotRGB(RGBstack,
        stretch = "lin")

# cir image
bands <- c(90, 34, 19)

CIRstack <- create_stack(f,
                         bands = bands,
                         epsg = epsg)

plotRGB(CIRstack,
        stretch = "lin")